#include <life/options.hpp>
#include <life/lifecore/application.hpp>
#include <life/lifealg/backend.hpp>
#include <life/lifediscr/functionspace.hpp>
#include <life/lifefilters/geotool.hpp>
#include <life/lifefilters/gmsh.hpp>
#include <life/lifefilters/exporter.hpp>
#include <life/lifevf/vf.hpp>


#include <iostream>

using namespace Life;
using namespace Life::vf;



namespace interp3d
{

    typedef Application Application_type;
    typedef boost::shared_ptr<Application_type> Application_ptrtype;

    enum type_champs
        {
            ScalarTest = 0,
            VectorialTest
        };

    inline
    po::options_description
    makeOptions()
    {
        po::options_description desc_options("test_interp_twomesh options");
        desc_options.add_options()
            ("hsize1", po::value<double>()->default_value( 0.1 ), "mesh size1")
            ("hsize2", po::value<double>()->default_value( 0.1 ), "mesh size2")
            ("gmsh", po::value<float>()->default_value( 2.1 ), " version of gmsh(2.0 or 2.1)")
            ;
        return desc_options.add( Life::life_options() );
    }

    /*_________________________________________________*
     * About
     *_________________________________________________*/

    inline
    AboutData
    makeAbout()
    {
        AboutData about( "Interp3d" ,
                         "Interp3d" ,
                         "0.1",
                         "Test1 verify geomap such as phi*ph^{-1}=Id ; Test2: on gamma, ||..v(u1)-..v(u2)|| < epsilon",
                         Life::AboutData::License_GPL,
                         "Copyright (c) 2010 Universite Joseph Fourier");

        about.addAuthor("Vincent Chabannes", "developer", "vincent.chabannes@imag.fr", "");
        return about;

    }


    /*_________________________________________________*
     * createMesh
     *_________________________________________________*/

    template<uint Dim,uint OrderGeo>
    boost::shared_ptr< Mesh<Simplex<Dim,OrderGeo,Dim> > >
    createMesh( std::string __name,std::string __str)
    {
        typedef Mesh<Simplex<Dim,OrderGeo,Dim> > mesh_type;
        typedef boost::shared_ptr<mesh_type> mesh_ptrtype;
        mesh_ptrtype mesh( new mesh_type );

        Gmsh gmsh;
        gmsh.setOrder(OrderGeo);

        std::string fname = gmsh.generate( __name, __str  );

        ImporterGmsh<mesh_type> import( fname );

        float gmsh_version = 2.1;//this->vm()["gmsh"].template as<float>();
        if (gmsh_version==2.0) import.setVersion( "2.0" );
        else if (gmsh_version==2.1) import.setVersion( "2.1" );

        mesh->accept( import );
        mesh->components().set ( MESH_RENUMBER|MESH_UPDATE_EDGES|MESH_UPDATE_FACES|MESH_CHECK );
        mesh->updateForUse();

        return mesh;

    }


    template<uint Dim,uint OrderGeo>
    boost::shared_ptr< Mesh<Simplex<Dim,OrderGeo,Dim> > >
    createMesh_hexaedre_p1( std::string __name,double meshSize)
    {

        GeoTool::Node x1(-1,-1,-1);
        GeoTool::Node x2( 1,-1,-1);
        GeoTool::Node x3( 1, 1,-1);
        GeoTool::Node x4(-1, 1,-1);
        GeoTool::Node x5(-1,-1,0);
        GeoTool::Node x6( 1,-1,0);
        GeoTool::Node x7( 1, 1,0);
        GeoTool::Node x8(-1, 1,0);
        GeoTool::Hexaedre H(meshSize,"Hexaedre",x1,x2,x3,x4,x5,x6,x7,x8);
        H.setMarker(_type="line",_name="Arete",_markerAll=true);
        H.setMarker(_type="surface",_name="Interface",_marker2=true);
        H.setMarker(_type="surface",_name="koko",
                    _marker5=true,
                    _marker3=true,
                    _marker4=true,
                    _marker6=true,
                    _marker1=true);


        H.setMarker(_type="volume",_name="Omega",_marker1=true);

        std::cout<<"\n" << H.geoStr() <<"\n";;
        typedef Mesh<Simplex<Dim,OrderGeo,Dim> > mesh_type;
        typedef boost::shared_ptr<mesh_type> mesh_ptrtype;
        mesh_ptrtype mesh( createMesh<Dim,OrderGeo>( __name + mesh_type::shape_type::name(),
                                                     H.geoStr()
                                                     /* +"Physical Volume(10) = {1};"*/ ));

        return mesh;

    }
    template<uint Dim,uint OrderGeo>
    boost::shared_ptr< Mesh<Simplex<Dim,OrderGeo,Dim> > >
    createMesh_hexaedre_p2( std::string __name,double meshSize)
    {
        GeoTool::Node x1(-1,-1,0);
        GeoTool::Node x2( 1,-1,0);
        GeoTool::Node x3( 1, 1,0);
        GeoTool::Node x4(-1, 1,0);
        GeoTool::Node x5(-1,-1,1);
        GeoTool::Node x6( 1,-1,1);
        GeoTool::Node x7( 1, 1,1);
        GeoTool::Node x8(-1, 1,1);
        GeoTool::Hexaedre H(meshSize,"Hexaedre",x1,x2,x3,x4,x5,x6,x7,x8);
        H.setMarker(_type="line",_name="Arete",_markerAll=true);
        H.setMarker(_type="surface",_name="Interface",_marker1=true);
        H.setMarker(_type="surface",_name="koko",
                    _marker2=true,
                    _marker5=true,
                    _marker4=true,
                    _marker3=true,
                    _marker6=true);


        H.setMarker(_type="volume",_name="Omega",_marker1=true);

        //std::cout<<"\n" << H.geoStr() <<"\n";;
        typedef Mesh<Simplex<Dim,OrderGeo,Dim> > mesh_type;
        typedef boost::shared_ptr<mesh_type> mesh_ptrtype;
        mesh_ptrtype mesh( createMesh<Dim,OrderGeo>( __name + mesh_type::shape_type::name(),
                                                     H.geoStr()
                                                     +"Physical Volume(10) = {1};" ));

        return mesh;

    }

    template<uint Dim,uint OrderGeo>
    boost::shared_ptr< Mesh<Simplex<Dim,OrderGeo,Dim> > >
    createMesh_rect_p1( std::string __name,double meshSize)
    {
        GeoTool::Node x1(-1,-1);
        GeoTool::Node x2( 1,-1);
        GeoTool::Node x3( 1, 1);
        GeoTool::Node x4(-1, 1);
        GeoTool::Quadrangle H(meshSize,"Rectangle",x1,x2,x3,x4);
        H.setMarker(_type="line",_name="Arete",_marker1=true,_marker3=true,_marker4=true);
        H.setMarker(_type="line",_name="Interface",_marker2=true);
        H.setMarker(_type="surface",_name="Omega",_marker1=true);

        //std::cout<<"\n" << H.geoStr() <<"\n";;
        typedef Mesh<Simplex<Dim,OrderGeo,Dim> > mesh_type;
        typedef boost::shared_ptr<mesh_type> mesh_ptrtype;
        mesh_ptrtype mesh( createMesh<Dim,OrderGeo>( __name + mesh_type::shape_type::name(),
                                                     H.geoStr() ) );

        return mesh;

    }

    template<typename element_type>
    void
    projectionH1(element_type & u,Application_ptrtype & test_app)
    {
        typedef Backend<double> backend_type;
        typedef boost::shared_ptr<backend_type> backend_ptrtype;
        typedef backend_type::sparse_matrix_type sparse_matrix_type;
        typedef backend_type::sparse_matrix_ptrtype sparse_matrix_ptrtype;
        typedef backend_type::vector_type vector_type;
        typedef backend_type::vector_ptrtype vector_ptrtype;

        typedef typename element_type::functionspace_type functionspace_type;
        typedef typename element_type::functionspace_ptrtype functionspace_ptrtype;

        static const uint16_type OrderChamp = functionspace_type::basis_type::nOrder;
        static const uint16_type OrderGeo = 1;

        std::cout <<"\ndefintion start\n";
        functionspace_ptrtype Xh=u.functionSpace();

        backend_ptrtype backend(backend_type::build( test_app->vm() ));

        sparse_matrix_ptrtype A = backend->newMatrix(Xh,Xh);
        vector_ptrtype F = backend->newVector(Xh);

        std::cout<<"\ndefinition finsih\n";

        auto pi = M_PI;
        //AUTO( Uex, cos(pi*Px()/2.)*cos(pi*Py()/2.) );
        //AUTO( grad_Uex , vec( -(pi/2.)*sin(pi*Px()/2.)*cos(pi*Py()/2.) , -(pi/2.)*cos(pi*Px()/2.)*sin(pi*Py()/2.) ) );
#if 1
        auto Uex = (Px()*Px()+Py()*Py())*(Pz()+2);
        auto grad_Uex = vec(2*Px()*Pz()+4*Px(),
                            2*Py()*Pz()+4*Py(),
                            Px()*Px()+Py()*Py() );
#else
        auto Uex = cos(pi*Px()/2.)*cos(pi*Py()/2.);
        auto grad_Uex = vec( -(pi/2.)*sin(pi*Px()/2.)*cos(pi*Py()/2.) , -(pi/2.)*cos(pi*Px()/2.)*sin(pi*Py()/2.) );
#endif

        std::cout<<"\nassemble1\n";
        form2( Xh, Xh, A, _init=true ) =
            integrate( elements(Xh->mesh()),_Q<2*(OrderChamp*OrderGeo)>(),
                                              idt(u)*id(u)
                                              +gradt(u)*trans(grad(u)) );
        std::cout<<"\nassemble2\n";

        form1( Xh, F, _init=true ) = integrate( elements(Xh->mesh()),
                                                _Q<15>(),//_Q<18>(),
                                                Uex*id(u)+grad(u)*grad_Uex );
        std::cout<<"\nassemble finsish\n";

        A->close();
        F->close();

        std::cout<<"\nsolve start\n";

        vector_ptrtype X= backend->newVector( Xh );
        backend->solve( A, A, X, F, false );
        u = *X;
        std::cout<<"\nsolve finsih\n";

    }

    template<uint Dim,uint OrderChamp,uint OrderGeo>
    void
    run(Application_ptrtype & test_app)
    {
        using namespace Life;
        using namespace Life::vf;


        Log() << "[testGeoMap] start test geomap inverse\n";

        //-----------------------------------------------------------------------------------//

        typedef Mesh<Simplex<Dim,OrderGeo,Dim> > mesh_type;
        typedef boost::shared_ptr<  mesh_type > mesh_ptrtype;

        typedef bases<Lagrange<OrderChamp,Scalar,Continuous,PointSetFekete> > basis_type;
        //typedef bases<Lagrange<OrderChamp,Scalar,PointSetFekete> > basis_type;
        typedef FunctionSpace<mesh_type, basis_type> space_type;
        typedef boost::shared_ptr<space_type> space_ptrtype;
        typedef typename space_type::element_type element_type;

        //-----------------------------------------------------------------------------------//

        std::cout<<"\nmesh_start\n";
        //double meshSize1=test_app->vm()["hsize1"].template as<double>();
        //double meshSize2=test_app->vm()["hsize2"].template as<double>();

            std::ostringstream __ostr;
            __ostr<<"graph_conv3d_P"<<OrderChamp<<".data";


        std::ofstream out;
        out.open( __ostr.str(), std::ios::out);

        double meshSize2=test_app->vm()["hsize2"].template as<double>();

        mesh_ptrtype mesh2 = createMesh_hexaedre_p1<Dim,OrderGeo>( "hexaedre2", meshSize2 );
        //mesh_ptrtype mesh2 = createMesh_rect_p1<Dim,OrderGeo>( "rectangle2", meshSize2);

        std::vector<double> __h(5);
        __h[0]=1;
        __h[1]=0.5;
        __h[2]=0.25;
        __h[3]=0.125;
        __h[4]=0.0625;

        typename std::vector<double>::const_iterator it = __h.begin();
        typename std::vector<double>::const_iterator it_end = __h.end();

        for( it; it!=it_end;++it)
            {
                double meshSize1 = *it;
                //double meshSize2 = *it;

                mesh_ptrtype mesh1 = createMesh_hexaedre_p1<Dim,OrderGeo>( "hexaedre1", meshSize1 );
                //mesh_ptrtype mesh1 = createMesh_rect_p1<Dim,OrderGeo>( "rectangle1", meshSize1);

                std::cout<<"\nspace_start\n";
                space_ptrtype Xh1 = space_type::New(mesh1);
                space_ptrtype Xh2 = space_type::New(mesh2);

                std::cout<<"\nelement_start\n";
                element_type u1(Xh1,"u1");
                element_type u2(Xh2,"u2");

                auto pi=M_PI;
                std::cout<<"\nprojection_start\n";
                //auto kk = cos(Px()*pi)*cos(Py()*pi)*cos(Pz()*pi);
                //auto kk = (Px()*Px()+Py()*Py())*(Pz()+2);
                //u1 = vf::project( Xh1, elements(mesh1), kk);
                //u2 = vf::project( Xh2, elements(mesh2), kk);
                //projectionH1<element_type>(u1,test_app);

                auto Uex = (Px()*Px()+Py()*Py())*(Pz()+2);
                auto grad_Uex = vec(2*Px()*Pz()+4*Px(),
                                    2*Py()*Pz()+4*Py(),
                                    Px()*Px()+Py()*Py() );

                projectionH1(u1,test_app);


                std::cout<<"\nintegration start\n";

                double errL2 = integrate( markedfaces(mesh2,mesh2->markerName("Interface") ),_Q<2*(OrderChamp*OrderGeo)>(),
                                                  (idv(u1)-Uex)*(idv(u1)-Uex) ).evaluate()(0,0);
                double errsemiH1 = integrate( markedfaces(mesh2,mesh2->markerName("Interface") ),_Q<2*(OrderChamp*OrderGeo)>(),
                                                    (gradv(u1)-trans(grad_Uex))*(trans(gradv(u1))-grad_Uex) ).evaluate()(0,0) ;

                std::cout << meshSize1 << " "<< std::sqrt(errL2) << " " << std::sqrt(errsemiH1) << "\n" ;
                out << meshSize1 << " "<< std::sqrt(errL2) << " " << std::sqrt(errsemiH1) << "\n" ;

            }


#if 0
                std::cout<<"\nexport_start\n";
                typedef Exporter<mesh_type> exporter_type;
                typedef boost::shared_ptr<exporter_type> exporter_ptrtype;
                exporter_ptrtype exporter( Exporter<mesh_type>::New( test_app->vm(), test_app->about().appName()+"1" ) );
                exporter->step(0)->setMesh( mesh1 );
                exporter->step(0)->add( "u1", u1 );
                exporter->save();
                typedef Exporter<mesh_type> exporter_type;
                typedef boost::shared_ptr<exporter_type> exporter_ptrtype;
                exporter_ptrtype exporter2( Exporter<mesh_type>::New( test_app->vm(), test_app->about().appName()+"2" ) );
                exporter2->step(0)->setMesh( mesh2 );
                exporter2->step(0)->add( "u2", u2 );
                exporter2->save();
#endif



    }

}// end namespace interp3d
int main(int argc, char** argv)
{

    interp3d::Application_ptrtype test_app(new interp3d::Application_type(argc, argv,
                                                      interp3d::makeAbout(),
                                                      interp3d::makeOptions()
                                                      ));

    test_app->changeRepository( boost::format( "/testsuite/lifediscr/%1%/" )
                                % test_app->about().appName()
                                );


    interp3d::run<3,1,1>(test_app);
    interp3d::run<3,2,1>(test_app);
    interp3d::run<3,3,1>(test_app);
    interp3d::run<3,4,1>(test_app);
    interp3d::run<3,5,1>(test_app);
    interp3d::run<3,6,1>(test_app);



    return 0;


}
